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We develop a new approach to building cosmological models, in which small pieces of perturbed 
Minkowski space are joined together at reflection-symmetric boundaries in order to form a global, 
dynamical space-time. Each piece of this patchwork universe is described using post-Newtonian grav¬ 
itational physics, with the large-scale expansion of the Universe being an emergent phenomenon. 
This approach to cosmology does not require any assumptions about non-local averaging processes. 
Our framework clarifies the relation between the weak-held limit of general relativity, and the cosmo¬ 
logical solutions that result from solving Einstein’s equations with a set of symmetry assumptions. 
It also allows the effects of structure formation on the large-scale expansion of the Universe to be in¬ 
vestigated without averaging anything. As an explicit example, we use this formalism to investigate 
the cosmological behaviour of a large number of regularly arranged point-like masses. In this case 
we hnd that the large-scale expansion is well modelled by a Eriedmann-like equation that contains 
terms that take the form of dust, radiation, and spatial curvature. The radiation term, while small 
compared to the dust term, is purely a result of the non-linearity of Einstein’s equations. 


I. INTRODUCTION 

The standard approach to cosmological modelling is a 
top-down one in which the first step is to solve for the 
homogeneous and isotropic large-scale expansion. Small 
fluctuations on large scales are then included using first- 
order perturbation theory [T], and large fluctuations on 
small scales are included by appealing to Newtonian the¬ 
ory m- This approach has many features that commend 
it as a good way to build cosmological models. Among 
the foremost of these is the mathematical simplicity in¬ 
volved at every step, as well as the fact that the resulting 
model has been found to be consistent with a wide array 
of cosmological observations (as long as dark matter and 
dark energy are allowed to be included). 

Nevertheless, while the top-down approach is simple 
and functional, it is not necessarily self-consistent or well- 
defined. This is because the standard approach assumes, 
from the outset, that the large-scale expansion of a statis¬ 
tically homogeneous and isotropic universe can be accu¬ 
rately modelled using a single homogeneous and isotropic 
solution of Einstein’s equations. It is far from obvious 
that such an assumption should be valid, as Einstein’s 
equations are non-linear, and because it has not yet been 
possible to find a unique and mathematically well-defined 
way of averaging tensors. This makes it extremely diffi¬ 
cult to assess the effect that the formation of structure 
has on the large-scale expansion of the Universe, with¬ 
out assuming that it is small from the outset. This is 
known as the “back-reaction” problem, which, despite 
much study, has uncertain consequences for the actual 
Universe [3]. 
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Part of the difficulty with the investigation of the back- 
reaction problem is that it is hard to do cosmology with¬ 
out first assuming that the geometry of the Universe can 
be treated (to at least a first approximation) as being a 
homogeneous and isotropic solution of Einstein’s equa¬ 
tions. This is, unfortunately, assuming the thing that 
one wants to question in the first place, which is ob¬ 
viously not ideal. The problems with studying back- 
reaction in the standard top-down approach to cosmo¬ 
logical modelling become increasingly apparent if one al¬ 
lows small-scale perturbations to the homogeneous and 
isotropic background. On small scales density contrasts 
must become highly non-linear, and extrapolation from 
the linear regime (which is assumed to be valid on large 
scales) can result in divergences [4]. On the other hand, 
appealing to the Newtonian theory results in a situation 
where the perturbations to the metric contribute terms to 
the field equations that are at least as large as the terms 
that come from the dynamical background, making the 
perturbative expansion itself poorly defined [5]. 

In this paper we report on an approach that side-steps 
all of these difficulties. We construct cosmological mod¬ 
els from the bottom up, by taking regions of perturbed 
Minkowski space, and patching them together using the 
appropriate junction conditions. We take the boundaries 
between these regions to be reflection symmetric, in or¬ 
der to make the problem tractable. The model that re¬ 
sults is a space-time that is periodic, and statistically 
homogeneous and isotropic on large scales, while being 
highly inhomogeneous and anisotropic on small scales. 
The equations that govern the large-scale expansion of 
the space are determined up to post-Newtonian accu¬ 
racy. The Newtonian-order equations reproduce the ex¬ 
pected behaviour of a Friedmann-Lemaitre-Robertson- 
Walker (FLRW) model filled with pressureless dust, while 
to post-Newtonian order we find new terms that appear 
in the effective Friedmann equations. For the case of a 
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single mass at the centre of each region, these terms all 
take the form of either dust or radiation. 


With a new generation of galaxy surveys on the hori¬ 
zon, that have the promise to map the structure in the 
Universe to unprecedented levels El m, it becomes in¬ 
creasingly important to understand the relativistic cor¬ 
rections that could be required in order to accurately 
interpret the data that results. To this end, higher-order 
corrections in perturbation theory are already being cal¬ 
culated (see, e.g., [8]). Our framework provides a way 
of consistently and simultaneously tracking the effects of 
relativistic gravity in both the regime of non-linear den¬ 
sity contrasts, and in the large-scale cosmological expan¬ 
sion. A proper understanding of these effects is required 
to ensure we understand all possible sources of error that 
could arise in the interpretation of the data, and also if 
we are to use the data to test and understand Einstein’s 
theory, and the dark components of the Universe. 

Our models do not rely on any assumptions about the 
large-scale expansion being well-modelled by any homo¬ 
geneous and isotropic solutions of Einstein’s equations. 
They are also well defined on small scales, as they are 
explicitly constructed from the post-Newtonian expan¬ 
sions that are routinely used to study the weak-field and 
slow-motion limit of general relativity. We are therefore 
able to model non-linear structure within the context of 
a cosmological model without falling foul of any of the 
problems outlined above (and that appear to be inherent 
in any top-down approach to cosmological modelling). 
The large-scale expansion of our model simply emerges, 
as a consequence of the junction conditions. We do not 
have to make any assumptions about the averaging of 
tensors, and do not have to assume anything about the 
existence of any background cosmology. 

Other recent approaches to bottom-up cosmological 
modelling include the application of geometrostatics to 
cosmology [9], as well as numerical relativity techniques 
pUi ITT] , perturbative approximation schemes [m El], 
and the re-discovery of the Lindquist-Wheeler models 
El- These studies have allowed the evolution of sub¬ 
spaces to be calculated [m ESI, numerical approxima¬ 
tions to both the space-time and the optical properties 
of the space-time to be determined USES], and the proof 
of interesting results such as the limit of many particles 
approaching a fluid [20] , and the non-perturbative nature 
of some structures [211 El] • Our study extends these pre¬ 
vious ones, we believe, by allowing increased flexibility 
for the distribution of matter, while maintaining a high 
degree of mathematical simplicity. 


In Section [n] we introduce the post-Newtonian formal¬ 
ism that we will use throughout the rest of the paper. In 
Section |m] we then explain how we will apply the junc¬ 
tion conditions, in order to build a cosmological model 
from regions described using the post-Newtonian approx¬ 
imation to gravity. Section [IV] contains a detailed presen¬ 
tation of the field equations and the junction conditions, 
both expanded to post-Newtonian levels of accuracy. In 
Section|V]we then manipulate these equations into a form 


that can be used to determine the motion of the bound¬ 
aries of each of our regions, and hence the expansion 
of the global space-time. We then proceed, in Section 
fVl| to explain how the general solution to the global ex¬ 
pansion can be calculated, for arbitrary distributions of 
matter, as well as for the special case of a single mass at 
the centre of each region. In both cases the lowest-order 
Newtonian-level solution to the equations of motion give 
a large-scale expansion that is similar to the dust domi¬ 
nated homogeneous and isotropic solutions of Einstein’s 
equations. To post-Newtonian order we find that, for 
the case of isolated masses, the only new terms that ap¬ 
pear in the effective Eriedmann equation look like dust 
and radiation. After that, in Section m we transform 
our solutions so that they are written as the evolution of 
proper distances in proper time, and on time-dependent 
backgrounds. We also consider the calculation of ob¬ 
servables in these models. Einally, in Section jvmj we 
conclude. 

We will use Greek letters (/i, z/, p, ...) to denote spatial 
indices, and Latin letters (a, 6, c, ...) to denote space- 
time indices. We reserve the latter part of the Latin al¬ 
phabet (i, jf, /c, ...) for the indices on the 2+I-dimensional 
hypersurfaces that will form the boundaries of each of 
our regions of space. Capital Latin letters (A, 5, C, ...) 
denote the spatial components of tensors on these bound¬ 
aries. As usual, a comma will be used to denote a partial 
derivative, such that 


where = t here is a time coordinate, and cp denotes any 
arbitrary function on space-time. Covariant derivatives 
will be represented by semi-colons. 


II. POST-NEWTONIAN FORMALISM 

The equations of General Relativity are known to re¬ 
duce to those of Newtonian gravity in the limit of slow 
motions {v <C c) and weak gravitational fields (^ ^ I). 
In the solar system, for example, gravity is weak enough 
for Newton’s theory to adequately explain almost all phe¬ 
nomena. However, there are certain effects that can only 
be explained using relativistic gravity. These include, for 
example, the shift in the perihelion of Mercury, which 
requires the use of relativistic gravity. To describe such 
situations it is useful to consider post-Newtonian gravi¬ 
tational physics. 

The post-Newtonian formalism is essentially based on 
small fluctuations around Minkowski space. Both the ge¬ 
ometry of space-time, and the components of the energy- 
momentum tensor, are then treated perturbatively, with 
an expansion parameter 


where t; = is the 3-velocity associated with the mat¬ 
ter fields, and c is the speed of light. The first step in 
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the post-Newtonian formalism is to associate all quanti¬ 
ties with an “order of smallness” in e. This is done for 
Newtonian and post-Newtonian gravitational potentials, 
as well as for every component of the energy-momentum 
tensor. 

In the remainder of this section we will outline the 
post-Newtonian expansion, and the quantities that are 
useful for solving the equations that result. Much of our 
discussion closely follows that of Will [23]. From this 
point on we will work in units where c = 1, unless explic¬ 
itly stated otherwise. 

A. Post-Newtonian Book-Keeping 


where v is the 3-velocity of the matter fields, and ip is any 
space and time dependent function in the system (such 
as one of the components of hab or Tab)- 

To find the lowest-order part of hu we note that the 
leading-order part of the equation of motion for a time¬ 
like particle takes the same form as in Newtonian theory. 
That is, vT I = ^ therefore have that 

the leading-order part of hu is 

hu ^ • (7) 

Similar considerations lead to the conclusion that the 
leading-order part of the spatial components of the met¬ 
ric are given by 


General relativity tells us that the curvature of space- 
time is directly related to its energy-momentum content 
by Einstein’s equations: 


hap ^ , ( 8 ) 

while those of the ta-components are given by 


Rab = SttG (^ab — -^Tga^ , (3) 

where Rab is the Ricci tensor, gab is the metric of 
space-time, G is Newton’s constant. Tab is the energy- 
momentum tensor (dependent on the matter content of 
the Universe), and T = g^^Tab is the trace of the energy- 
momentum tensor. These are a set of ten non-linear par¬ 
tial differential equations in four variables. 

In the vicinity of weakly gravitating systems we take 
the metric to be given by 


9ab — Vab T hab ■> (4) 

where r]ab = diag{ — l^ 1) is the metric of Minkowski 

space, and hab are perturbations to that metric. We also 
take the energy-momentum tensor to be given by 

(5) 


where g is the energy density of the matter fields mea¬ 
sured by an observer following p is the isotropic 
pressure, and is a time-like unit 4-vector, given by 
where r is the proper time along the inte¬ 
gral curves of and where is normalized such that 
u^Ua = —1. Anisotropic pressure could have been in¬ 
cluded in Eq. , but would only appear at O(e^) in 
gap, and O(e^) in the equation of motion of time-like 
particles. This is beyond the level of accuracy used in 
this paper, and so we do not include it here. We can now 
expand hab, hi P and in orders of e, and relate the 
resultant quantities to each other via Eq. (§. 

To begin this we first note that, in the post-Newtonian 
formalism, time derivatives add an extra degree of small¬ 
ness to the object they operate on, as compared to spatial 
derivatives. This follows because the time variations of 
the metric and energy-momentum tensors are taken to 
be a result of the motion of the matter in the space-time, 
such that 


hta - e" ■ (9) 

The next-to-leading-order parts of each of these compo¬ 
nents is 0{e^) smaller than the leading-order part, in 
every case. 

Similarly, the lowest-order part of g can be determined 
from the leading-order part of the tt-component of Eq. 
(©■ This takes the form of the Newton-Poisson equation, 
= —SirGg^ so that the lowest-order part of g can 
be seen to be 


g rsj . (10) 

Here, and throughout, we have chosen units such that 
spatial derivatives do not change the order-of-smallness 
of the object on which they operate. To find the lowest- 
order contribution to the pressure we can consider the 
conservation of energy-momentum, = 0. The 

lowest-order part of the spatial component of these equa¬ 
tions is g{u^)^t + huh(u^) p = ^ghu,a —p,ai from which 
it can be seen that 


p ^ . (11) 

Again, the next-to-leading-order part of the energy den¬ 
sity is O(e^) smaller than the leading-order part, while 
the higher-order corrections to the pressure will not be 
required for what follows. 


B. Post-Newtonian Potentials 

In order to solve the equations of both Newtonian and 
post-Newtonian gravitational physics it is useful to define 
some potentials, as well as make some identifications for 
the components of the energy-momentum tensor. The 
first of these involves the leading-order part of the energy 
density, which we write as 


(6) 




( 2 ) 


= P, 


<P,t ^ </^,7 > 


( 12 ) 
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where p is the density of mass. Here, and throughout, 
a superscript in brackets denotes a quantity’s order-of- 
smallness in e. The next-to-leading order part of the 
energy density is then written as 

= (13) 

where H is known as the the specific energy density. In 
what follows, we will also use to denote the spatial 
components of at lowest order. 

Using p we can define the first of our potentials, which 
is simply the Newtonian gravitational potential, defined 
implicitly as the solution to 

=-47rG/), (14) 

where = dada is the 3-dimensional Laplacian oper- 
ator of flat space. The reason for defining the potential 
in this way is very simple: it allows us to write the solu¬ 
tion to h^t ,77 ~ —SttGp as = 24>. At this point the 
“solution” for hu is little more than a change of notation 
(with obvious historical significance). When it comes to 
post-Newtonian potentials, however, the equations be¬ 
come much more complicated. This change of notation 
is then much more useful, especially if the potentials that 
we define are simply solutions to Poisson’s equations. 

With this in mind, it is useful to make the following 
implicit definitions for new gravitational potentials 

V\ = -2^, 

= -AnGpv^ , 

V^4>i = —AirGpv^ , 

V^4>2 = —47rGp4>, 

V^4>3 = -47rGpn, 

V^4>4 = —AirGp , 

where v‘^ = y ^ e^, ^ and 4>i ^ 4>2 ^ 

$4 ~ e^- In what follows, we will not require any 
potentials of order higher than e^. 


C. Green’s Functions 

When applying the post-Newtonian formalism to iso¬ 
lated gravitational systems, such as the Sun, it is usual to 
assume that the space-time is asymptotically flat. This 
allows the boundary terms of Green’s functions to be ne¬ 
glected, so that the solutions to Poisson’s equation are 
given by simple integrals over spatial volumes. This sim¬ 
plicity is a substantial benefit when solving for the po¬ 
tentials defined in the equations above. In the case of a 
cosmological model, however, we cannot assume asymp¬ 
totic flatness. We must, therefore, be more careful with 
the boundary terms. 

The equations we want to solve take the form of the 
Poisson equation, 

VV = -^, (15) 


where if is now being used to denote a potential, and 
where T is some function on space-time (such as the mass 
density). To solve this equation we consider the Green’s 
function, ^(x, x',t), that satisfies 




-(5(x - x') + Gi, 


(16) 


where x and x' denote spatial positions, where S(x — x') 
is the Dirac delta function, and where Gi is a constant 
over spatial hypersurfaces (the need for Gi in this equa¬ 
tion will become apparent shortly). 

We want to solve Eq. © over a spatial volume U, 
with boundary dfl. The Green’s function we will use 
for this must, of course, satisfy Gauss’ theorem on this 
domain, such that 


f dV = [ n-VGdS, (17) 

Jn JdQ 

where n is the unit-vector normal to the boundary. If we 
now choose n • = 0 as the boundary condition for 

Q, then Eqs. ( [l^ and © imply 




( 18 ) 


where V is the spatial volume of the cell. The solution 
to Eq. (15) is then given, in terms of 0, by considering 


f QT dV = [ gV^<fdV 

= / [V ■ (gv^) - wg ■ v</j] dv 

Jq 

= / [V ■ (gv</?) - V ■ (<fivg) + dv 

Jn 

= / V ' {QVp - pVQ) dV - ip(p, 

Jn 

where (p = Gi f^cp dV = y dV is a constant over 
U. Rearranging, the potential (p can be seen to be given 
by 


p = p- 


[ gj^dV-i- [ Gn-VpdA, (19) 

Jn JdQ 

where we have again made use of the boundary condition 

^ = 0 . 

If one were now to assume that p was asymptotically 
flat, then the first and last terms on the right-hand side 
of Eq. (19) would vanish. The Green’s function would 
then take the form of the Newton kernel, such that 


If fp 

(^(x7) = -— / I--dV. 

47r |x-x'| 


( 20 ) 


However, these assumptions are not true in general, and 
especially not in the case of cosmological modelling. This 
is because cosmological models do not have asymptoti¬ 
cally flat regions, by their very definition. In what fol¬ 
lows, we must therefore be more careful. We have to 
specify appropriate boundary conditions for our poten¬ 
tials, and include the boundary terms in Eq. (19), if we 


are to use the Green’s function formalism to determine 
the form of our gravitational potentials. 
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III. BUILDING A COSMOLOGY USING 
JUNCTION CONDITIONS 

In this paper we will take a bottom-up approach to cos¬ 
mological modelling. This will involve considering cosmo¬ 
logical models that are constructed from large numbers 
of cells, that can be put next to each other to form a 
periodic lattice structure. The shape of each cell will be 
taken to be a regular polyhedron, and will be assumed 
to be identical to every other cell, up to rotations, reflec¬ 
tions and translations. 

The physical systems that we intend to model with 
these cells will depend on the size of cell that we are con¬ 
sidering. For example, for cells that are approximately 
the size of the homogeneity scale (about 100 Mpc), we 
could consider modelling clusters of galaxies, as illus¬ 
trated in Fig. [2 Other systems, such as individual 
galaxies, could equally well be modelled with cell sizes 
of the order of about 1 Mpc. The only requirement we 
have is that the system must satisfy the requirements of 
the post-Newtonian formalism. Specifically, this means 
that V c and p ^ p, so that the bulk of the interior of 
each cell is described by Newtonian and post-Newtonian 
gravitational physics. 

The post-Newtonian formalism is expected to work 
well in the regime of non-linear density contrasts, and 
so should be expected to be adequate for modelling most 
aspects of the gravitational fields of galaxies and clusters 
of galaxies. This formalism is, however, limited to scales 
much smaller than the cosmological horizon. We there¬ 
fore require each of our cells to be much smaller than 
the Hubble radius, A violation of this requirement 

would result in matter at the boundary of a cell moving 
at close to the speed of light. We also assume each cell 
is filled with normal matter, so that pressures are small 
with respect to energy densities. 

We note that the post-Newtonian framework cannot, 
and should not, be used to describe multiple cells simul¬ 
taneously. However, due to the periodicity of our lattice 
structure, we only need to know the space-time geometry 
of any one cell, and its boundary conditions with neigh¬ 
bouring cells. As we will show below, this information 
is sufficient to tell us how we should expect the entire 
Universe to evolve. 

Let us now turn to a more detailed consideration of 
the conditions required in order to join two cells together 
at a boundary. First and foremost, the cells must satisfy 
certain smoothness requirements across their respective 
boundaries, known as the Israel junction conditions, if 
their union is to be a solution to Einstein’s equations. 
These conditions, in the absence of surface layers, are 
given by [24] 

bij] = 0 (21) 

[Kij]=0, ( 22 ) 

where [ip] = — (p^~^ denotes the jump across the 

boundary for any quantity (p, and the i and j indices de¬ 
note tensor components on the boundary. The and 



FIG. 1. Two adjacent cubic cells, with example matter 
content, consisting of filaments and voids. The second cell is 
the mirror image of the first. This figure was produced using 
an image from [25] . 


superscripts here show that a quantity is to be eval¬ 
uated on either side of the boundary (i.e. on the sides 
labelled by + or —, respectively). 

In these equations, is the induced metric on the 
boundary, and Kij is the extrinsic curvature of the 
boundary, defined by 


_ dx^ 


and 


_ dx^ dx^ 


(23) 


(24) 


where denotes the coordinates on the boundary, and 
is the space-like unit vector normal to the boundary. 

In our construction we choose to consider reflection 
symmetric boundaries. The Israel junction conditions 
can then be simplified. The situation we wish to con¬ 
sider is illustrated in Fig. [^ for two cubic cells. We 
use x^ and x^ to denote the coordinates used within the 
first and second cells, respectively. Reflection symmetry 
means that Eq. (21) is automatically satisfied. The sec¬ 
ond junction condition, given by Eq. (22), can be written 
as 


dx^ dx^ dx^ dx^ (_) , 

where and n~ ^ are outward and inward pointing 
normals, in the first and second cells, respectively. They 
are shown in Eig. Now, mirror symmetry implies that 
Symmetry therefore demands that 

dx^ dx^ (+) dx^ dx^ (+) 

This implies that Kij = —Kij, or, in other words, that 
the extrinsic curvature must vanish on the boundary of 
every cell, i.e. 


Kij=0. 


(27) 
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FIG. 2. A schematic diagram showing the normal vectors 
involved in the junction conditions, and . The vector 
is shown as a dashed arrow, and is the mirror image of 


This equation must be satisfied on each and every bound¬ 
ary in our lattice. 

Eq. (27) is valid for any cell shape, as long as the 
boundaries are reflection symmetric. In general, how¬ 
ever, there are only a finite number of ways that spaces of 
constant curvature can be tiled with regular convex poly- 
hedra. These are listed in Table [ij where we have also 
included the Schlafii symbols that allow the structure of 
the lattice to be inferred, and the number of cells in each 
of the different structures. A Schlafii symbol {pqr} cor¬ 
responds to a lattice with p edges to every cell face, q cell 
faces meeting at every vertex of every cell, and r cells 
meeting around every cell edge. One may note that reg¬ 
ular lattices on hyper-spherical spaces have a maximum 
of 600 cells, while those on flat and hyperbolic lattices 
always have an infinite number of cells. 

In what follows we will often require knowledge of the 
total surface area of a cell, A. For each of the five differ¬ 
ent polyhedra in Table [l| we can write A = where 


Cell 

Faces per 

Surface Area 

Shape 

Cell, K 

Coefficients, 

Tetrahedron 

4 

24^3 

Cube 

6 

24 

Octahedron 

8 

12^3 

Dodecahedron 

Icosahedron 

12 

20 

-,on V 25 +IOV 5 
25+11^5 

120\/3 

7+3v^ 


TABLE II. The five possible different cell shapes, together 
with the number of faces per cell, and the numerical coeffi¬ 
cients for the surface area, a,^ = AjX^. 


X is the distance from the centre of a cell to the centre of 
a cell face, and where k denotes the number of faces per 
cell (e.g. a cube has k = 6 faces). For the purposes of 
creating the diagram in Fig. we chose to consider cube¬ 
shaped cells. In what follows we will also often consider 
this particular case. One of the advantages of tessellat- 
ing the Universe into cubic cells is that tilings of this 
type exist for open, closed and fiat universes. Another 
advantage is that Cartesian coordinates can be used, and 
aligned with the symmetries of the cells. 

Let us now finish this section by briefly considering 
the motion of a boundary at x = X{t^y,z)^ where the 
x-direction has been chosen to be orthogonal to the cen¬ 
tre of a cell face (as in Fig. [^. The 4-velocity of this 
boundary is then given by the following expression 


Ua 


dx^ 

dr 


dt 

dr 



(28) 


Lattice 

Structure 

Lattice 

Curvature 

Cell 

Shape 

Cells per 

Lattice 

{333} 

+ 

Tetrahedron 

5 

{433} 

+ 

Cube 

8 

{334} 

+ 

Tetrahedron 

16 

{343} 

+ 

Octahedron 

24 

{533} 

+ 

Dodecahedron 

120 

{335} 

+ 

Tetrahedron 

600 

{434} 

0 

Cube 

00 

{435} 

- 

Cube 

00 

{534} 

- 

Dodecahedron 

00 

{535} 

- 

Dodecahedron 

00 

{353} 

- 

Icosahedron 

00 


where x^ are the coordinates of points on the boundary, 
where r is the proper time measured along a time-like 
curve in the boundary, and where we have chosen the 
integral curves of to stay at fixed y and ^ coordinates. 
This vector is orthogonal to the space-like normal to the 
cell face, such that U^Ua = 0. 

We can now define two types of derivatives along the 
boundary, in time-like and space-like directions, respec¬ 
tively. These are given to lowest order by 

■ = u^da = dt^x^td, 

1 ^ = m'^da = dt + X^aOx , (29) 

where is a space-like vector in the cell face, which 
satisfies m^ria = 0. These expressions allow us to write 
the lowest-order parts of rit and ua as 


TABLE 1. A summary of all regular lattice structures that 
can exist on 3-surfaces of constant curvature. Hyper-spherical 
lattices are denoted by +, flat lattices by 0, and hyperbolic 
lattices by —. The lattice structure is given by the Schlafii 
symbols, {pqr}, which are explained in the text. The shape of 
the cells, and the number of cells in the lattice are also given. 
For further details of these structures see [26] . 


nt = -rixX^t 

nA = -rixX^A ■ (30) 

We know that Xt ^ e and ^ 1, which means that 
the leading-order part of ^ e. This information will 
be used below. 
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IV. GOVERNING EQUATIONS and 


In this section we will present the equations that gov¬ 
ern the dynamics of each cell, and its matter content. We 
will first use Einstein’s field equations to relate the space- 
time metric (given in Eq. @) to the energy-momentum 
content of the cell (given in Eq. ^). After this, we 
will evaluate the extrinsic curvature of the cell bound¬ 
aries, using the same geometry. The extrinsic curvature 
will be required to vanish, in order to satisfy the reflec¬ 
tion symmetric boundary conditions, and will provide us 
with the information required to study the evolution of 
the boundary. 


A. Einstein’s Field Equations 


To begin this study, we need to evaluate the Ricci ten¬ 
sor for the perturbed Minkowski space given in Eq. 
Recall that the leading-order contributions to the metric 
perturbations are at hu ^ ^ e^, and ht^ ^ e^. The 

leading-order components of the Ricci tensor are then 
given by 



Here we have used the short-hand notation = 9^9^, 

and before, we have chosen 

units such that spatial derivatives do not add an order 
of smallness, and have used super-scripts in brackets to 
denote the order of a quantity. 

The only higher-order part of the Ricci tensor that we 
require will be the 0{e^) part of the tt-component, which 
is given by 


d(4) _ £ 




jiv/.: 


( 2 ) ,2 

tt \ 






(34) 




No other components of the Ricci tensor will be required 
at this order. 

To proceed further we now need to make a gauge 
choice, in order to eliminate superfluous degrees of free¬ 
dom. Eor this we use the standard post-Newtonian 
gauge, which is given by [23] 


^(3) _i^(2) _Q 


(36) 


Note that this is not the same as the harmonic gauge. 

The perturbed metric, Eq. Q, and the definition of 
the 4-velocity can now be used to write the 4-velocity as 


(l + ^ + -)(l;t,A‘)+0(e4) 


(37) 


where = v^v^. We can use this equation to give us 
the components of up to post-Newtonian levels of 
accuracy. We note that in order to evaluate the field 
equations at O(e^) we only need to know Tu = —T = p. 
Using Eqs. and ( j^ , and the gauge conditions ( j^ 
and (36), the field equations § then give us 

= —87rGp + 0(e^) , (38) 

+ O(e^) . (39) 

Using the potentials defined in Section |n| we then find 

(40) 

and 




2$(5 


flU • 


(41) 


These solutions, together with our gauge conditions 


(35) and (36), allow us to simplify Eq. (|33|) , to get 
A3) _ 


R 


'tfj, 


/7^3) ^ 


The field equations ^ then give 


= iQ-n-GpVf,, 


which has the solution 


'h/i — 


-4U^ + -X,t ^, 


(42) 


(43) 


(44) 


where we have again made use of the potentials defined in 
Sectionim Eqs. (40), (41) and (44) give the leading-order 
contributions to all of the components of the perturbed 
metric. 

To go further, we now need to evaluate . This will 
be done using Eq. (34), our gauge conditions (35) and 
(36), and the lowest-order solutions found above. The 
relevant part of the Ricci tensor then simplifies to 


d(4) _ 

rCf-i- — 


1 


- 4|V$p - 


(45) 


which, using the identity |V4>y = — 4>V^4>, can 

be written as 


h(2) , ,(2) _ q(2) _n 


(35) 


--- 


V^(2$^) -8i>V^i> +V^/i; 


2^(4) 


(46) 























Similarly, we can write the tt-component of the right- 
hand side of Eq. © as 


lTgtt=p{v^-^ + ln+^P^]+0{e^) . (47) 


Equating Eqs. (46) and (47), and using the field equa¬ 
tions ([^, we then find that 


— —24*^ -|- 44*1 H“ 44*2 H“ 24*3 64*4 . 


(48) 


Once more, this solution has been written in terms of the 
potentials defined in Section [HJ Eqs. (40), ( [41] ), (44) and 
(48) give all of the components of the metric that we will 


require. 


B. Extrinsic Curvature Equations 


Let us now calculate the extrinsic curvature of the 
boundary. To do this we require the covariant deriva¬ 
tive of the normal, ria-iy. As stated earlier, the leading 
order parts of rit and are 0(e) and 0(1), respectively. 
The components of can then be seen to be given, up 

to the required order, by 




n 


Il'flU fHt,u 








(49) 


+ O(e^), 


and 




T 


(50) 




- (3) 


+ 0(e"): 


and 




(51) 


“ + 0{e^) ■ 


New lines have been used, in each of these equations, to 
separate terms of different orders. 

The extrinsic curvature of the cell boundaries can now 
be calculated using Eq. (24). The tt-component of this 
equation is given, to lowest order, by 




n 


(52) 


where we have used rit = —UxX^f At next-to-leading- 


order have 

V‘2 ^(2) 

Kit> = ^n.(eu - 


-2X 


htt^x^xX^t 


_L _ h(->> I ') 

I 2 'hx,v I 'Hv,x) 

(3) 


n^x'’X,tt 
(3) , .(3) 




fht.v 


2 ' 


(53) 


These equations can be simplified even further by making 
use of the result ua = —UxX^a- 

Similarly, the leading-order parts of tA and AB- 
components of the extrinsic curvature tensor are given 
by 

, (54) 

and 

k21 = -X,ab ■ (55) 

These two equations, together with the result Kij = 0, 
imply that X^a is independent of t, ^ and ^ at lowest or¬ 
der. This implies that X^a also vanishes at lowest order, 
as X^A is forced by symmetry to vanish at the centre of 
every cell face. 

This information allows us to write simplified versions 
of the next-to-leading-order parts of the tA and AB- 
components of the extrinsic curvature tensor as 


A- - 

^tA — ^,At ^ 9 


r(2) 


'Ha,x 'hx,A 


h^2\Xt 

(56) 


and 


+ 7;nahAB,a • 


(57) 


This is all the information we require about the extrinsic 
curvature of the boundaries in our lattices. 

Using Eq. (27), we can finally obtain from Eqs. (52) - 

.(4) 


(57) the conditions 


X,u = 


and 


hi 






(58) 

+o(ey, 




= o 


hilx-ht]A-^^,AX,t 


0{e^), (59) 


cc=X 


and 


X^AB = ^AB^,x 


x=X 


+ Oie^), 


(60) 


where we have made explicit the requirement that each 
equation is to be evaluated on the boundary, at x = X. 
Eq. ( [ 5 ^ is similar to the geodesic equation, as shown in 
Appendix 
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V. COSMOLOGICAL EXPANSION 

We now have enough information to find the equations 
for the acceleration of the boundary of each cell, up to 
post-Newtonian accuracy. These will be the analogue of 
the Friedmann equations, of homogeneous and isotropic 
cosmological models. In this section, we begin by re¬ 
producing the lowest order Friedmann-like equations at 
Newtonian order. We then proceed to obtain the post- 
Newtonian contributions to the same equations. 


where C = C{y^ z) is an integration constant in t. How¬ 
ever, if this equation is to satisfy X,tA = O(e^), we see 
that C must also be a constant in y and Eq. (66) 


is similar in form to the Friedmann equation, with X 
behaving like the scale factor, and the constant C behav¬ 
ing like the Gaussian curvature of homogeneous spatial 
sections. 

From now on, we will take the positive branch in Eq. 


(66), which corresponds to an expanding universe. The 


solution to this equation depends on the sign of C. Eor 
C = 0 we have 


A. Newtonian Order 


We can begin by defining the gravitational mass within 
each cell as 

M = [ p , (61) 

Jn 

where is the spatial volume element at zeroth order, 

and Q is now the spatial volume of the interior of a cell. 
We can apply Gauss’ theorem to this equation, so that 


47rGM = - 


Jn 

[ n • V4> dA^^^ , 

JdQ 


(62) 


X = 



/ ISttGM 


t-to 


2/3 


0(e2 


(67) 


where to is a constant, which can be absorbed into t by 
a coordinate re-definition. 

Eor G 7 ^ 0, we can obtain parametric solutions. Eor 
G > 0 we have 


X = 


SttGM 

a.\C\ 


sin^ 



t-to 


4ttGM , . , 


( 68 ) 


where n = ria is the normal to a cell face ,and dA is the 
area element of the boundary of the cell, dCl. 

By noting that the cell face is flat to lowest order (i.e. 
that X A = 0(e^)), we can see that it is possible to write 
X = X(t). Together with the lowest-order part of Eq. 
( [^ , this implies that n • V4> is constant on the boundary. 
We therefore have 


where y = f dt/X is the analogue of conformal time. 
Similarly, for G < 0, we get 


X = 


SttGM 


sinh^ 



AttGM = —An • V4>, 


(63) 


where A is the total surface area of a cell. Eor a cell that 
is a regular polyhedron, we can take the surface area to 
be A = where X is the coordinate distance from 

the centre of the cell to the centre of the cell face, where 
K is the number of faces of the cell, and where are the 
numerical coefficients that are given in Table [HI 


The lowest-order part of Eq. (58), along with Eq. (63), 
then gives us 


X,, = (n-V^)U=x = - 


47rGM 

A 


AttGM 

' C^kX2 


(64) 


We can solve this equation by multiplying both sides by 
X^f This gives 


l((Y 4 :TTGMX^t 

which can be integrated to find 


(65) 


X, 


± 


SttGM 


C, 


( 66 ) 


t-to 


AttGM 


(69) 


These solutions represent parabolic, closed, and hyper¬ 
bolic spaces, respectively. At this order, they expand 
in the same way as a dust-dominated homogeneous- 
and-isotropic geometry with the same total gravitational 
mass. 


B. Post-Newtonian Order 


To find the post-Newtonian contributions to the accel¬ 
eration of the boundary of each cell it will be useful to 
replace all of the terms in Eq. (58) by more physically 
relevant quantities. The first step in doing this will be to 
include all orders in the calculation of Gauss’ theorem, up 
to O(e^). We can begin by expanding the normal deriva¬ 
tive of the potential on the boundary to post-Newtonian 
order, so that 
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1 

2 


n -Vh+t dA = 




c(2) 






+ K [ + 0(e®), 

Js 


(70) 


where dS is the area element of one face of the cell, where 
hz is the number of faces of the cell, where S is the area 
of one face of the cell, where A = kS is the total surface 
area of the cell, and where dS^‘^^ is the O(e^) part of the 
area element. 

We can then apply Gauss’ theorem, making sure we 
use covariant derivatives, in order to ensure we include 
all post-Newtonian terms. This gives 

i / n • Vhtt dA 

= \ ( dV 

^ Jn 

= J dV + 0(e®). (71) 

( 2 ) 

Note that here hu is being used to denote both hi^ and 
Expanding, and using the lower-order parts of the 
field equations, allows us to write the first and last terms 
on the right-hand side of Eq. as 

- [ V^htt dV = -AttG f p dV^°'> - AttG [ p dV^^'> 

2 Js Jn Jn 

+ dV^°^+0{e^), (72) 

and 

J \W^\^dV = J (^47rGi>/)+ dV^°'> 

= 4ttG{p^) + k [ dS^°'>, (73) 

Js 

where dV^‘^^ is the O(e^) correction to the volume ele¬ 
ment, and where we have introduced the new notation 


{ip)= [ ifdV, (74) 

Jn 

where if is some scalar function on the space-time. 

To find the area and volume elements up to O(e^) we 

[ 2 ] 

need the induced 2-metric on the boundary, , and the 
spatial 3-metric, both up to O(e^). These are given 


by 

[2] [3] 

9ab = 9ab - 

= (1 + 2<^)Sab + O(e^) , (75) 

and 

= (1 + ‘2^)S^u + O(e^). (76) 

The determinants of these two quantities, up to the re¬ 
quired accuracy, are given by 

d^t{9^AB) = 1 + 44> + 0{e ^), (77) 

and 

det{9]^l) = 1 + 64> -f O(e ^). (78) 

By taking the square root of these determinants, and 
Taylor expanding them, we obtain the higher-order area 
and volume elements in terms of their lower-order coun¬ 
terparts: 


= 2$ dS^°'>, (79) 

^ 3^ ^y(0) _ (gQ) 


We can now evaluate the higher-order corrections to the 
normal. 

As X^A vanishes at lowest order, is given by 

^a(2) ^ _^(2) _ 

We can also use the normalisation of the space-like nor¬ 
mal, rio^rJ^ = 1, to obtain 


,^(2) _ 


= - 4 >. 


(82) 


Using Eqs. (70)-(82) we can now write 


Ht,x 


^ / 

Js 

= - AttGM +l [ +2 k [ 

^ Jn Js ' 

+ K [ dS^°K (83) 

Js ’ 


To understand this equation further, we note that the 
second term on the right-hand side can be written as 
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- [ = [ (-V2$2 + 2V2$i + 2V2$2 + + 3V^^4) dV^°'> 

2 Jn Jn 

= - [ n- ^^( 0 ) f ( 2 v 2 $i + 2 V 2$2 + + 3V2$4) dV^°'> 

Js Jn 

= -2k [ dS^°'> + [ (2V2$i + 2V2$2 + V2$3 + 3V2$4) dV^°'> , 

Js ’ 


(84) 


r 


where we have now used Eq. (48) and Gauss’ theorem. where we have used Eqs. (85) and (88), and substituted 


Using Eqs. (74) and (84), we can now re-write Eq. (83) 
as 

.(4) 


K, 


L 




h 


'tt,X 




dS 


= — AttGM — S7rG{pv‘^) — 87rG(/)4>) 
- 47rG(pn) - UirGip) • 


(85) 


To proceed further, let us now determine the functional 
form of X up to O(e^). Using the lowest order parts of 
Eqs. (58) and (pM, this is given by 


^ = at) + \{y^ + a)n ■ + 0{eS , (86) 


where ({t) is some function of time only, 
derivatives, and substituting from Eq. (|^| 


Taking time 
, then gives 


C« =X,tt - liy^ + z^in ■ Vi>)- + 0(e®) 


=$ r - 




^(4) 

_ l(3) _ 
2 'hx,t 


Xt - 34>.tX, 


4^4 


-{y^ + za{n-W^)-+0{a), (87) 


where all terms in this equation should be taken as being 
evaluated on the boundary. 

Many of the terms in this equation can be simplified 
using the lower-order solutions. Eor example, using Eqs. 
and ( [6^ , and taking time derivatives, gives 

, 22A7t‘^G‘^M‘^ 247rGMG , , 

(«■ v$)- =-• (88) 


a 


a, 




We are now in a position to express the equation of mo¬ 
tion in terms of variables that can be easily associated 
with the matter fields in the space-time. Recall that the 
total surface area of the cell at lowest order is given by 
A = kS = a«(X(0))2, j^(o) is the zeroth-order 

part of X, which we solved for earlier. As Eq. (87) is 


a function of t only, we can integrate over the area on a 
cell face to obtain 


^C,tt 

= - AttGM 


ii 


a^Xz \ a^X J 

( 87rGM$ 


~ titx,t ~ 1 dS 


— 8'!tG{pv'^) — 87rG(p$) — 47rG(/9n) — 127rG(p) 
nUn'^G'^M'^ 12 ttGMC\ 




dS 


in for lower-order solutions. 


We can make use of our gauge condition, hty^y = 
\hyy^t = and Gauss’ theorem, to replace one of 

the terms in this equation in the following way: 


■L 


^ I — 3 / ^.tt dV . 


> / ^,t 

Jn 


(90) 


Using Eqs. (90), we can write Eq. (89) as 


^C,tt 

= - AnGM + 


kS 


■L 


a«(X(0))2 
87rGM$ 


96tt^G^M^ 

a,.X(o) 


- 12nGMC j 


- 


( 0 ) 


dS-3 [ dV 

Jq 


,S V««(X(0)) 

— 8'kG{pv^) — 8'kG{p^) — inGlpTl) — 127rG(p) 
/1127r2G2M2 127rGMG\ 


”V al{X(0))5 a^{X(0))4 


)Ia 


-h z^) dS + O(e^). 
(91) 


Then we must find the O(e^) correction to the surface 
area. This is to ensure that we include all 0{€^) cor¬ 
rections to the equation of motion. To find the O(e^) 
correction to the surface area, we must take into account 
that the edges of the boundaries are curved at O(e^) and 
the background is also curved at O(e^). Then the surface 
area is given by 


A = sS = K, 




2$) dS + O(e^) 


(92) 


In general, this is dependent on the cell shape. Therefore, 
the equation of motion of the boundaries can be written 
in its final form as 


f 0(e®): 
(89) 
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X 


AttGM UttGMC 


AttG 


,tt 


+ 


+ 


A a,(X(0))2 a^(X(0))2 L' 

K f f SttGM^ 


2{pv‘^) + 2(pi>) + (pH) + 3{p) 


I 




/ 

Jq 


dV + 




a,(X(o))2 Js W{X(o))^ 

flUn^G'^M^ 127rGMC\r k f.^ 2 . ,a ^ 2 2^1 6^ 

( a2(.Y<«))'. - aJxmAlaJXfolf JJ’' +Oj+0(0^ 


gGTT^G^M^ 


(93) 


This gives us the post-Newtonian correction to the accel¬ 
eration of the boundary, and hence the post-Newtonian 
correction to the accelerated expansion of the Universe. 
This equation is the main result of this section. 

The first term in Eq. (93) is the standard Friedmann- 
like term for a dust-like source. The second term is a 
higher-order correction due to the presence of the spatial 
curvature-like term. The third term in the first line con¬ 
tains all post-Newtonian corrections to the matter sec¬ 
tor. In the second line, the terms integrated over the 
area and volume are dependent on the potential, and 
the rate of change of the potential. The terms that go as 
X~^ behave like radiation terms, when compared to the 
standard Friedmann equation. However, we remind the 
reader that these terms are purely a result of geometry, 
and the non-linearity of Einstein’s equations. In the last 
line of this equation, the post-Newtonian contributions 
depend on the cell shape that is being considered. 


VI. POST-NEWTONIAN COSMOLOGICAL 
SOLUTIONS 


In this section we will present solutions to Eq. (93). 


We start by specialising this equation to lattices con¬ 
structed from cubic cells. In this case we have = 24^ 
and K = 6. The total surface area is given by 


A = 24C - 


WttGMX 


i 


+ 12 / ^dS + O(e^) (94) 


where ( is the time-dependent part of X. The lowest- 
order part of the limits of integration in both the y and 
^-directions are also given simply by —and X^^\ 
The equation of motion, given in Eq. (93), then simplifies 
to 


X 


ttG 


M + 5MG + ‘2.{pv^) + 2{p^) + (pH) + 3(p) 


+ 


8 ( X ( 0))2 ( 


+ 


(X(0))3 




27 


fTw^G^M'^ ttGMG 
V36(X(0))5 “ 2(X(o))4 


f A^'kGM 
3(X(0))2 


1 dS • 


/ 

Jn 


dV 


(y^ + z^) + 0{e^). 


(95) 


The last term in this equation is a function of its position 
on the boundary, and vanishes at the centre of a cell face. 

We can solve Eq. ( [9^ if we know the functional form of 
the potential, and its time dependence, as well as that 
of the post-Newtonian corrections to the matter sector. 
However, we do not need to know the functional form 
of all the post-Newtonian corrections, as we can replace 
one of the higher-order terms with the conserved post- 
Newtonian mass. This is given by 


Mpn = p(^^v^ + 3$^ dV = ^{pv^) + 3(p$). (96) 

The proof that this object is conserved can be found in 
Appendix 

We can now find the general functional form of the 
potential ^ for our model, using the Green’s function 
formalism. We will do this below for the case of cubic 
cells. Similar analyses can also be performed for the other 
platonic solids. This result can then be used to evaluate 
the acceleration of the boundary. After this, we proceed 
to study the special case of point sources, where the form 
of the potential ^ can be found somewhat more straight¬ 
forwardly, and where the first post-Newtonian correction 
to the acceleration of the boundary can be determined 
explicitly. 


A. The General Solution: An Application of the 
Green’s Function Formalism 


In this section we wi ll use the explanation of Green’s 

and in particular Eq. (19). 


functions from Section HC 


To give a concrete example of how this works, let us now 
consider a lattice of cubes arranged on In Figure 

we show a 2-D representation of such a 3-D lattice. 
As before, we assume reflection symmetry about every 
boundary, which imposes a periodicity on our structure. 
For cubic cells of edge length L = 2X, the periodicity of 
the the lattice will be 2L. That is, if we move a distance 
of 2L in any direction in our lattice, then two reflections 
will ensure we return to a point that is identical to our 
starting position. 

If we consider one example cell, then we can now use 
the method of images to construct a Green’s function 
that is symmetric around each of its boundaries, and 
that therefore satisfies the required boundary condition 
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2L 




m 

# 

• 

2L 

t 



f 

Y 

\x'W 


and its images. The source for the Green’s function is 
then the sum of the Dirac combs that contain all 8 of the 
image points described by these 8 position vectors, plus 
an additive time dependent constant. Using Eqs. (16) 
and (97), the source function is thus given by 


i=i 


1 


p7ri/3- 




1 


(98) 


E E 8£,3' 


• ^ (x-x'(j)) 

— z —^ 


i=i /3ez2 


(99) 


where does not include the null triplet /3 = (0,0,0), 
and where the are given by 


= x' 


= x' + Lei — 2^1 ei , 


FIG. 3. A 2-D representation of the vectors used to locate 
the position of image points. In two dimensions we require 
only four unique vectors, as compared to eight in 3-D. The 
four lattice vectors are given by = x', x'^^^ = x' + Lei — 
2x'iei, x'^^^ = x'-hLe 2 — 2 x 262 , and x'^^^ = —x'-hL(ei + 62 ). 


= x' + Le2 — 2^262 , 

= -x' + L(ei + 62) + 2x363 , 


at each of its faces, n-= 0. Due to the identical na¬ 
ture of every cell, such a Green’s function can then be re¬ 
used for each of the cells. The way that this method will 
work is by introducing mirror images of the points in our 
original cell. We therefore consider the point sources of 
the Green’s function to be a set of Dirac delta functions, 
separated from infinitely many identical point sources by 
pairwise distances of 2L. The structure that results will 
be a superposition of several ‘Dirac combs’. 

A Dirac comb can be expressed as a Fourier series in 
the following way: 

(97) 

/3ez3 /3ez3 

where (3 = (/3i, / 32 , /ds), and where p 2 and are inte¬ 
gers. To construct our Green’s function, we must include 
the location of image points, in relation to the location of 
points in the central cell. In Figure in 2-D, we choose 
an arbitrary point in the example cell, and use x*^^^ to 
represent its position with respect to the centre of that 
cell. The mirror symmetry across the boundary of every 
cell results in 8 image points in the 8 surrounding cells. 
However, we only require 4 unique vectors to describe the 
positions of the initial point source and its images. These 
4 vectors are shown in Figure The other points can be 
added by considering Dirac combs, with periodicity 2 L, 
that contain these initial 4 points. 

In 3-D we can do something similar, but we require 8 
unique vectors to describe the position of the initial point 


= x' + Le 3 - 2 x 363 , 

= -x' + L( 6 i + 63 ) + 2 x 362 , 
= -x' + L (62 + 63 ) + 2 Xi 6 i , 


x'^«) = -x' - 


■ ^(^1 + 62 + 63 ) , 


( 100 ) 


where x'^ = x' • 61, X2 = x' • 62 and X 3 = x' • 63, and 
where 61, 62 and 63 are orthogonal unit vectors. 

The solution to Eq. (99) is then given by 




1 


j=i /3ez3 


87r^L|/3|^ 


^7Tif3- 




( 101 ) 


If we take, as an example equation, the Newton-Poisson 
equation (14), we can then see from Eq. (19) that the 


Newtonian potential 4> is given by 


^ = ^ + 47rG [ GpdV - 

Jn 


ttGM 




( 102 ) 


where we have used Eq. (64). This is the general so¬ 


lution for the potential. That is, for any given energy 
density distribution we can simply evaluate the integral 
to find the potential, as well as the rate of change of the 
potential. In Eq. (95), this, along with the mass, velocity 


and pressure of the matter fields can be used to evaluate 
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the acceleration of the boundary numerically, up to post- 
Newtonian accuracy. Hence, this expression allows us to 
find the post-Newtonian correction to the expansion of 
the Universe for general configurations of matter. 

Of course, within the interior of each cell we can also 
solve for each of the post-Newtonian potentials defined 
in Section |IIB| using the same Green’s function, as each 
of these potentials is defined as the solution to a Pois¬ 
son equation. We therefore have a complete solution, for 
the motion of the boundary of every cell, and for the 
geometry interior to each cell. 


B. A Special Case: Point Sources 


In order to find an even more explicit solution, let us 
now consider the case of a single point mass, located at 
the centre of each cell. In this case, the Poisson equation 
simplifies to 


= -4nGM5{x) , (103) 


where M is the mass we defined in Eq. We can 

again use the method of images to solve for 4>. In this 
case we want to place our image points such that 4> satis¬ 
fies the inhomogeneous boundary condition given in Eq. 
(64). We therefore place an image mass at the centre 
of each surrounding cell, so that image masses are sep¬ 
arated by a distance L from each other. We can then 
express the source of the potential, 4>, as a sum of Dirac 
delta functions that correspond to these masses. We then 
continue by placing image masses at the centre of every 
cell that surrounds the cells that already contain image 
masses (taking care not to place two masses in any one 
given cell). We repeat this process Af times, and then let 
Af ^ oo. 


This description may initially sound similar to the pro¬ 
cess used to find the Green’s function, above. There is, 
however, a subtle difference. In order for there to be a 
non-zero normal derivative of on the boundary, we need 
to take a sum whose number of terms tends to infinity, 
rather than an array that is infinitely extended from the 
outset. These two things are not equivalent, in this case. 

The source of the potential can then be written as 


J\f 

= -AttGM lim V 5{x - 1/3) , (104) 


where L = 2X is again the edge length of the cubic cell, 
and w here A/" is a positive integer. The solution to Eq. 
(104) is given by 


4> = lim 

A/'^c 


N 

E 

' 0=-M 


GM 

|x-i/ 3 | 


fit) > 


(105) 


where /(t) is an arbitrary function of time. 

We can use f{t) to regularize the value of 4> at x = 0, 
such that it reduces to the regular form for a Newtonian 
potential around a single point source. This can be done 
by subtracting the contribution of all image points to the 
potential at x = 0: 


AT 


4> = 


lim 

Af^oo 


E 

I3=-M 


GM 

|x-i/ 3 | 


M 


lim 

Af^oo 



GM 


(106) 


where /3* indicates that the null triplet f5 = (0,0,0) is 
excluded from the sum. With this choice, the value of 
4> near x = 0 does not change as the number of image 
masses is increased. This can be considered as a bound¬ 
ary condition imposed at the location of the mass. _ 

As L is the only time dependent quantity in Eq. (106), 
it can be seen that the rate of change of 4> is simply given 

by 


Af^oo 

(3* = -Af 


(107) 


+ lim / 

Af^oo ^^ ^ 
f3=-Af 


3GM{f3 • 


|x-/ 3 L |5 


GML 2GM. 

- hm > , ^ ’ -hm > ——— 

^ ^ ISIL^ Af^oo ^ ^ 1311 

/3*=-Ar i3*=-Ar 


Now, in order to solve Eq. ( [9^ , we need to evaluate 
a few integrals. We need 4> and 4> ^ integrated over the 
boundary of a cell, and integrated over the volume. 
These integrals are given explicitly below. Eirstly, 



AT 


lim 

Af^oo 


E 

I3=-N 


GM 

y/{L/2 - + (y - P 2 LY + {z- 


Af 


lim 

Af^oo 


E 

l3-=-M 


GM' 

m. 


dydz . 
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This can be simplified by redefining coordinates such that y = Ly, and z = Lz. This gives 

/ LI2 pl/2 r A/" 

$U=L /2 dydz = GML / lim V 

-i/2 i-1/2 » 

= GMDL , 


AT 


_1_ 

' V (1/2 - /3i)2 + (y - /32)2 + (5 - /33)2 l/3|. 


dydz 


(108) 


where the last line of this equation defines the quantity D. Secondly, integrated over the boundary is given by 

rL/2 f>Ll2 r n-h/r( ^ ^ _\Pl\2t^t ^ 


/ iv/z _ 

^.tU=L /2 dydz = / lim W 

-L/2 J-L/2 I - , 


GMiP^x + p2y + p3Z-\f3\^L)L,t 


+ lim 

\f —^nn ^ ^ 


GML^t 


l-L/2 J-L/2 [(^ “ PlLY + (^ - hLy + (^ - |/3|T2 

We can again simplify this integral by using y and z coordinates. This gives 

(^i(l/2) + I 32 y + - |/ 3 p) 


dydz . 


nL/2 .1/2 

/ ^Ax=l /2 dydz = GML^t 

J-L/2 J- 


M 

lim 

Kf —VPiT) ‘ ^ 


N ^ 

lim > T— 

Kf —Vnn ^ ^ /S 


[( 1/2 - + (// - / 32)2 + (5 - / 33 ) 2]^/2 \(d\ 


dydz 

(109) 


/- 1/2 

= GMEL^t , 

where the last line gives the definition of E. Finally, the second time derivative of <1>, integrated over the volume, is 

.L/2 

/ dxdydz 

J-L/2 

— r P2y E f^z — \l3f LYGM[pix^21/E fsz — \fdf L)L^tt 

J-L/2W^^^^j^ |x-/3L|5 |x-/3L|3 


- E 

|x-/9Ep 


, GMitt , 2GML/^ 

lim > ’ - lim > ’ 




=gmL, 


/• 1/2 

J-1/2 


1 - 1/2 L 

GMLL 


lim V 3(/3ix + /? 2 / + feg-|/3p)^ _ ^ 

f —Vrv^ Y ^ A/* —Vrv^ ^ ^ 


^*=-M 
N 


dxdydz 




i3=-Ar 


r E 

E, 


|X-^|5 
{I3ix + l32y + Psz - |/3|2) 


l/^P _ ]• w — 

A/’^oo ^^ lx — /3|^ A/’^oo ^ 1/31 

I3=-Ar ' /3* = -A/^ 1^1-' 


dxdydz 


AT 


■ lim T 

i^iJ 


0=-N 


|x-^|= 


dxdydz 


=GML^F + GMLL^ttP , 

where F and P are defined by the last line. 


For a single point source at the centre of a cell w e take 
= p = n = (p^) = 0. Making use of Eqs. (108), (109) 
and (110), as well as the lower-order solutions, we then 


(110) 


I- 

volving and to find 

GM 


find that Eq. (95) reduces to 


X 


,tt 


X,tt = - 


GM 


QC^ 


TT + 57rC - 9EC - 3FC 


6X2 
ttG^M^ 
6X3 


TT + 57rC - 9EG - 3FC 


2D + - - F -3E 


IdTT 

“ 9 “ 


2D+--F-3E 


6(X(o))3 

ttGMC 


IdTT 

~9' 




+ 0{e^). 

( 112 ) 


V36(X(o))5 2 (X(o))4 


(/ + ^2) + 0(e®) 


( 111 ) 


To solve this equation, we can evaluate it at the centre 
of the cell face, where y = 0 and z = 0. In this case, 

( = X and it is convenient to recombine the terms in- 


where the last line defines N and J. 

As can be seen from Table |ml and the plots in Fig¬ 
ure]^ the numerical constants D^E^E and P are all of 
order unity, and converge rapidly as the number of im¬ 
age masses becom es la rge. The two terms on the right- 
hand side of Eq. (112) look like dust and radiation, re¬ 
spectively. However, J = 1.277rG^M^ is positive, so the 
radiation-like term appears to have negative energy den- 















































16 




FIG. 4. The percentage difference from the asymptotic value 
of D, E, F and P, for various numbers of image points in the 
partial sum. 


FIG. 5. Illustrations of the solutions of X(t) for different val¬ 
ues of C and J. The constant J will likely take much smaller 
values, for realistic configurations. Its value is exaggerated 
here to illustrate the effect of the post-Newtonian terms. 


sity, and hence contributes positively to the acceleration 
of the boundary. _ 

Integrating Eq. (112) gives the Friedmann-like equa¬ 
tion 


Xl = ^-E_c + 0{e^), (113) 

where C is a constant. The solutions to Eq. ( |113[ ) depend 
on the value of C. If (7 = 0 then 




J 

m ’ 


If 0 < C < ^ then 

X=^± JCsin [VCti] , 

G G 


be positive for there to be real solutions, which is always 

TV t2 

true if J > 0 and C < 0. Eor C ^ ^ there are no real 
solutions. In Eigurej^we present the functional form of 
X (t) for three example values of C and J. 

If we were to extrapolate these solutions beyond their 
reasonable regime of validity, to very early times, then we 
would observe a bounce at the following minimum values 
of X: 




N VN'^-JC 
C C ’ 

J 

2N^ 


if C 7^0 
if C = 0 . 


(117) 


This concludes our discussion of explicit solutions for 
X(t), in the presence of point-like particles. 


VII. RELATIONSHIP WITH FLRW MODELS, 
AND OBSERVABLES 


t-to = - JCcos [VCrj] . (115) 

In both of these equations 77 = J dt/X is analogous to the 
conformal time parameter used in Eriedmann cosmology. 
Eor C <0, on the other hand, we obtain 


± (t - to) 
N 


In 


- X + CX + V- 


(-C)3/2 

V-CX2 ^ 2XX - J 


+ 


(-C) 


-CX2 + 2XX - J 
(116) 


where we have written the inverted function t(X), for 
convenience. The arguments under the square root must 


So far, we have only calculated coordinate distances, in 
terms of coordinate time. In this section we will trans¬ 
form these quantities into proper distances and proper 


Constant 

Asymptotic Value 

D 

1.44... 

E 

0.643... 

F 

-1.63... 

P 

0.304... 


TABLE III. The asymptotic values of P, P, F and P, which 
are approached as the number of image masses diverges to 
inhnity. 
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time. We will then perform a coordinate transforma¬ 
tion that allows the space-time within each cell to be 
written as perturbations on a homogeneous and isotropic 
Robertson-Walker background. Finally, we will work out 
the expansion rates of our cell edges in this new de¬ 
scription. This will allow us the clearest possible com¬ 
parison with standard Friedmann-Lemaitre-Robertson- 
Walker (FLRW) cosmological models. As an example, 
we will sometimes use the model described in Section 


VIB Similar analyses can be performed for other config¬ 


urations, by adjusting the calculations that follow. 


A. The Proper Length of Cell Edges 

For our example cubic cell, we choose to consider an 
edge defined by the intersection of cell faces at x = 
X(t^y^z) and y = Y{t,x^z). The proper length of such 
a curve, in a hypersurface of constant t, is given by 

£ = / y/(l -1- 2 ^){dx^ + dy‘^ + dz‘^) + O(e^). ( 118 ) 

J edge 

Expanding this expression, and taking the locations of 
the relevant corners of the cell to be at 2 ; = ±Zc, then 
gives 


£ = 


[ \l + ^) dz + 0{e^). (119) 

J-Zc 


Using Eqs. (63) and (86), it can be seen that 

ttGM 


^ 2 


6 


+ 0(6^), 


( 120 ) 


where L is being used here to denote the coordinate dis¬ 
tance between the centres of two cell faces, on opposite 
sides of our cubic cell. 


Using the solution derived in Section VIB 


we can nu¬ 


merically evaluated the in tegra l in Eq. (119), between 
the limits specified in Eq. (120), to find 


£-L-0.125GM. 


( 121 ) 


Using this result, the Eriedmann-like equation (113) can 
then be written as 


where all terms are to be evaluated at the corner of the 
cell, and where the observer has been taken to be moving 
with velocity equal to in each of the x, y and z- 

directions. This equation can be re-written using the 
solution in Section [VI B[ to find 

. , 3.61GM 3G 

^ - 1 +-^-^ . (124) 


where we have used Eq. (121) to replace coordinate dis¬ 


tanc es wi th proper length. In terms of this proper time, 
Eq. ( 122[ ) is given by 


dC 

dr 


(16V - 54.0GMG) 4.A1G‘^M‘^ 


C 


C? 


- 4G(1 - 3G). 
(125) 


The funct ional form of this equation is again the same as 
Eq. (113), but now with corrections to both the dust-like 


and spatial curvature-like terms, as well as the radiation¬ 
like term. 


B. Transformation to a Time-dependent 
Background 


Up to this point we have treated the geometry within 
each cell as a perturbation to Minkowski space. This 
is convenient, as it allows the apparatus constructed for 
the study of post-Newtonian gravity in isolated systems 
to be applied with only minimal modifications. When 
considering the case of a cosmological model, however, it 
is also of use to be able to understand the form of the 
gravitational fields in the background of a Robertson- 
Walker geometry. In this section we will show that 
the perturbed Minkowski space description and the per¬ 
turbed Robertson-Walker description are isometric to 
each other, as long as we restrict our consideration to 
a single cell. We will present coordinate transformations 
that make this isometry explicit. This discussion follows, 
and extends, that presented in m- 

We begin by writing the unperturbed line-element for 
a Robertson-Walker geometry as 


ds‘^ = —dP + a(t)^ 


(ddP + dy‘^ + dP) 

[l + |(^2+^2+^2)]2 ’ 


(126) 


- 

— 


16 V 


64.9G^M‘^ 

£2 


-4G, 


( 122 ) 


The functional form of this equation is obviously unal¬ 
tered, with only a small additional contribution to the 
radiation-like term. 

Let us now write this equation in terms of proper time, 
r, of an observer following the boundary at the corner 
of the cell. In this case coordinate time can be related 
to proper time using the normalisation of the 4-velocity 
tangent to the boundary, U^Ua = —1, so that 




dt 

dr 


= l + $ + 5(Xy + 0(e4), 


where a{i) and k are the scale factor and curvature of hy¬ 
persurfaces of constant t, respectively. Erom the Eried- 
mann equations we know that k ^ (<^ ^ ^ it' 

a ^ 1, this means that we must immediately require that 
/c ^ e^, which means that k can be treated as a (homo¬ 
geneous) perturbation to a spatially-flat background. 

With this information, we can write the line-element 
of a perturbed Robertson-Walker geometry as 

ds^ — (1 — — h^^^)dP + 2a{t)lP?^dx^di 

+ a(tf + U?l - +y'^ + 5 ^)^ dx^^dx '', 

( 127 ) 


(123) 
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where the hatted coordinates = (x^y^z), and where 
are perturbations. In this expression we have writ¬ 
ten the tt-component of the metric to O(e^), the flu- 
component to 0(6^), and the t/i-component to O(e^). 


The lowest-order part of the Einstei n’s f ield equations 
give, for the perturbed metric in Eq. (127), 


h 






(128) 


and 

SirGpa^ + = —6a = 3(a £)^a + 3ka . (129) 


An integrability condition of this latter equation is that 
SnGpa^ + = constant = 2 C 4 , (130) 


where by constant we mean not a function of t, x, y or 


z. This result shows that the term in Eq. (|129) 


behaves like dust, in the way that it sources the evolution 
of the scale factor. 

Let us now consider the following coordinate transfor¬ 
mations: 




X 

X = — 

a 




+ 0(e") 


y = 




0(6-) 


a 


1 + + 


+ 0(e4), (131) 


where T{t^x^) is an unspecified quantity of O(e^). The 
scale factor on the right-hand side of these equations is 
written as a function of the time coordinate t, and is 
related to a{t) by 


a{i) = a(t) 




+ 0{e^). (132) 


In the im-hatted coordinates, the perturbed metric in 
Eq. (127) transforms into the one given in Eq. 
where the perturbations around the Robertson-Walker 
geometry are given in terms of the perturbations about 


Minkowski space in the following way: 

hff = hlf - + y^ + Z^) + 0(6^^) 


hi"!} +T^ + 2^x^h 


^tfJ, 


^tfJ, 




^ ^ ^2 ^ ^2) ^ 0(6^) 


2C4 _k\ a-, 

^ ~ 2/a3' 

= 4 ") + 2Tt + 2'^{hfl + T^)x^^ 

/5C'4^ kCi U\, 2 2 2 x 2 6x 

(133) 

These equations cannot be used to transform a global 
perturbed Robertson-Walker geometry to a global per¬ 
turbed Minkowski geometry, as the velocities that result 
in the latter would be greater than the speed of light on 
scales of . They are, however, perfectly sufficient 
to transform the geometry within any one of our cells. 
As every cell is identical, this provides us with a way to 
transform the entire geometry. 

We can use these transformations to relate the proper 
length calculated in the Minkowski space background, £, 
to the proper length in a flat Robertson-Walker back¬ 
ground, C. The relevant corners in this background are 
given hy z = ±Zc, which can be related to the positions 
of the cell corn ers in the Minkowski space background 
using Eq. (131), so that 


^c = — + 

a 


4a^ 


+ 0(e") 


(134) 


Eor a flat Robertson-Walker background we do not nee d 
any O(e^) corrections to Zc, as we did for Zc in Eq. (120), 
because the boundaries are flat to this order m- 


Using Eqs. (131) and (133), C is given by 


C'Zi J a{ 1 ^ ) dz 

y* a ^1 - (x^ -h y^ z^) ^ ] dz . (135) 


(^ 4 )^ ,2 , ^ 2 \ 


Using ^ = a^ at lowest order, and again choosing the 
edge of our cell nt x = X and y = Y, we can integrate 
part of this equation to find 


C- 


1: 









3 


3 



4> dz . 


(136) 

The last term in this equation can be evaluated numeri¬ 
cally, in the same way we evaluated the last term in Eq. 
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(121). We will do this below, for the explicit solution 


found in Section |VIB[ consisting of a single point-like 
mass at the centre of every cell. 

Using Eqs. (134) and (|136|) we then obtain 


C-L- 


~£- 


ttGM 5X^{a^tf 
~3 
5X^C4 
18a3 


12a2 


0.922GM 


(137) 


where we have used Eqs. (121) and (129) in the last 


line. For a flat FLRW background we can use Eq. (129), 


together with the solution for the scale factor, to find 




(138) 


where C 4 has been taken to be positive. Now, using Eqs. 
(67) and (138), for the solutions of X and a, respectively. 


we can infer that 


C4 = 


nGMa^ 

2X3 


(139) 


The proper length of a cell edge in the Robertson-Walker 
background is therefore given by 


Cc^C- 0.436GM . 


(140) 


The Eriedmann-like equation can then be written in 
terms of this quantity as 



(leX - bd.OGMC) SMG^M^ 

I ^2 


- 4C(1 - 3C), 
(141) 


This equation has the same form as Eq. (125), but with a 


different numerical coefficient for the radiation-like term. 
Again, we remind the reader that this term, although 
it takes the form of a radiation fluid in the evolution 
equation for the scale of the space, does not correspond 
to any actual matter field. It is purely a result of the 
non-linearity of Einstein’s equations. 


C. Observables 

The majority of this paper has been aimed at con¬ 
structing a geometry to model the space-time of a uni¬ 
verse that can accommodate large non-linear density con¬ 
trasts on small scales, while maintaining a periodic struc¬ 
ture that means it is invariant under a certain set of dis¬ 
crete spatial translations on large scales. Space-time ge¬ 
ometry itself, however, is not a direct observable. We 
will therefore pause, and consider in this section the con¬ 
sequences of this new type of model for cosmological ob¬ 
servables. 

We will briefly consider two types of observations: (i) 
those based on the motion of time-like particles (i.e. 


galaxies, dark matter, etc.), and (ii) those based on ob¬ 
servations made along null geodesics. Of course, almost 
all observations are in fact made by observing light, but 
one can in principle split ones considerations into the ef¬ 
fect of the expansion of space on the light directly, and 
the effect of the expansion on the astrophysical objects 
that emit or interact with the light. 

Let us start with the motion of time-like particles. The 
model constructed above, and in particular the Green’s 
function formalism adopted in Section |VI A[ allows n- 
body simulations of structure formation to be modelled 
to post-Newtonian accuracy in a fully self-consistent way. 
This problem has recently been studied using a differ¬ 
ent approach in m, where actual simulations were per¬ 
formed, and observables extracted. The effects of post- 
Newtonian gravitational fields are found to be small in 
all aspects of these simulations, but are nonetheless im¬ 
portant for understanding the consequences of relativis¬ 
tic gravity, and in order to accurately interpret precision 
cosmological data. We expect our formalism to produce 
similar results, although the simplicity of our equations 
may bring some computational benefits. 

Winding back the evolution of the Universe to earlier 
times, the consequences of relativistic gravitational fields 
should be expected to become more significant. The cor¬ 
rections we found for the large-scale expansion of our 
particular solution, presented in Section ra were of 
the order of our cell size squared divided by the Hubble 
radius squared. Taking cells the size of the homogeneity 
scale, this gives corrections at the level of about 1 part in 
10^. Larger cells will give larger corrections. By coinci¬ 
dence, this is about the same size as the contribution of 
Cosmic Microwave Background (CMB) radiation to the 
expansion of the Universe at late-times. Our corrections 
also scale like a radiation fluid, so naively extrapolating 
our results suggests that our corrections may start to be¬ 
come significant in the evolution of the Universe at about 
the same time as radiation starts to become important. 
For high precision observables, such as the CMB, it is 
conceivable that this could have some impact on the in¬ 
terpretation of data. 

Let us now consider observations based on the evolu¬ 
tion of null congruences in the space-time. The method 
of propagating null geodesics through a lattice universe 
of this type has already been considered in [13]. The 
method suggested there was to decompose the tangent 
vector to these curves as 


+ rf'). 


(142) 


where and are an orthogonal pair of time-like and 
space-like unit vectors. When a null geodesic reaches a 
cell boundary one can then read off u^ka and n^, which 
can be transformed to similar quantities in the new cell 
using the coordinate transformations from Eqs. (23) and 


These quantities can then be used as initial con¬ 
ditions for extending the same null geodesic through the 
new cell. 

Depending on the coordinate system being used, the 
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cosmological redshift in these models comes primarily 
either from the expansion of space or the motion of 
the cell boundaries (for the time-dependent and time- 
independent backgrounds, respectively). As the leading 
order part of the expansion of each of our cells is given 
by a Friedmann-like equation with a dust-like content, we 
expect this to contribute the usual leading-order term to 
the redshifts calculated over large distances. What our 
framework adds beyond the usual treatment, however, is 
a consistent way to calculate the consequences of both 
Newtonian and post-Newtonian gravitational potentials 
on all scales up to the size of the cells being considered 
(larger scales can always be included using standard cos¬ 
mological perturbation theory). The effect of relativistic 
gravity on galaxy redshifts has recently attracted much 
attention, in the light of upcoming surveys [8]. 

Finally, the effect of inhomogeneity on the brightness 
and shear of a source can be calculated using the Sachs 
optical equations [28] : 

^ + a*a = 

^ + 2a0 = Cabcd'm°‘k'’m'^k'‘‘ 

dX 

^ + 2w6> = 0, (143) 

cLX 

where a and uj are the expansion, shear, and vorticity 
scalars of the rays of light that are being modelled. The 
affine distance along these rays is denoted A, and Rah 
and Cahcd are the Ricci and Weyl tensors that describe 
the curvature of the space-time. The complex vectors 
are a set of mutually orthogonal and normalised vectors 
on the screen space orthogonal to and the angular 
diameter distance can be calculated by integrating the 
expansion scalar, such that va ^ exp(f^ OdA). 

One can immediately notice that, when light propa¬ 
gates in the near-vacuum regions between galaxies, the 
driving term in the evolution equation for d is absent, 
while the driving term in the corresponding equation for 
a is non-zero. This is exactly the opposite of what hap¬ 
pens in homogeneous and isotropic space-times [29] . Our 
model provides a way in which this switching on and off of 
terms can be consistently modelled to post-Newtonian or¬ 
der in an expanding universe. Such quantities are of vital 
importance not only for constructing Hubble diagrams, 
but also for correctly interpreting data from galaxy sur¬ 
veys and the CMB. It is therefore important that the 
models used to interpret these observations are as robust 
as possible. 


VIII. CONCLUSIONS 


We have constructed cosmological models by patching 
together many sub-horizon-sized regions of space, each 
of which is described using post-Newtonian physics. The 
boundaries between each of these regions was assumed 
to be reflection symmetric, in order to make the problem 
tractable. This allowed us to find the general form for 
the equation of motion of the boundary, as well as the 
general form of the post-Newtonian gravitational fields 
that arise for general matter content. These results both 
follow from a straightforward application of the junction 
conditions, which in the latter case provides the appropri¬ 
ate boundary conditions for solving Einstein’s field equa¬ 
tions. 

As an example of how to apply this formalism, we con¬ 
sidered the case of a large number of isolated masses, each 
of which is positioned at the centre of a cubic cell. We 
found that the large-scale evolution that emerges from 
such a configuration is well modelled by an equation that 
looks very much like the Friedmann equation of gen¬ 
eral relativity, with pressureless dust and radiation as 
sources. The radiation term is necessarily much smaller 
than the dust term (if the post-Newtonian expansion is 
to be valid), and appears in the Friedmann-like equation 
as if it had a negative energy density. This happens with¬ 
out any violation of the energy conditions, as the term in 
question arises from the non-linearity of Einstein’s equa¬ 
tions, and does not directly correspond to any matter 
content. 

While small, the radiation-like term appears as the first 
relativistic correction to the large-scale expansion of the 
Universe, when the matter content is arranged in the 
way described above. This term provides a small nega¬ 
tive contribution to the rate of expansion, and a small 
positive contribution to the rate of acceleration. The ex¬ 
istence of a radiation-like term has been found previously 
using exact results derived for the evolution of reflection 
symmetric boundaries m, and using the shortwave ap¬ 
proximation for fluctuations around a background metric 

m- 


It remains to be seen what form the first relativistic 
correction will take for more general configurations of 
matter fields, or for structures built with less restrictive 
symmetries. The former of these cases can be studied 
using the approach we have prescribed in Section |VI A[ 
while the latter requires a generalisation of our formal¬ 
ism. 
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Appendix A: Geodesic Equation 


Appendix B: Post-Newtonian Mass 


The geodesic equation, in terms of proper time along 
a curve, is given by 


(Px^ dx^ dx^ 


where T^^^ are the Christoffel symbols. Let us now con¬ 
sider the motion of a boundary at x = X{t^y^z). The 
proper time derivatives of this boundary are given by 


and 


d^X 

dr‘^ 


dr ’ 




(A2) 


+ X ABX^X^r + 


,rr • 


(A3) 


Eq. (A2) then allows us to write the t, y and 2 : compo¬ 


nents of the geodesic equation in terms of partial deriva¬ 
tives as 

t,TT = + 0(e®) , (A4) 

x,TT = + 0(e®) , (A5) 

■ (A6) 

Using these equations, Eq. ( |A3| ) can be written as 
(fX 


dT‘^ 


= {t,rY 


X,u - XA\eAtX"t + 


+ X^ABX^tX^t - A 




x=X 

(A7) 


The Christoffel symbols required for evaluating Eq. ( |A1| ) 
can be simplified using Eq. (A7), and written explicitly 
as 


= + 0(e"), 


(AS) 

(A9) 


YX b c 


^ a, + 2<1>^ ; 


J- 1, 

— + k 


( 3 ) 

tx,t 


+ 2i>,tA ( + 2^^^X/ + 0(e®), (AlO) 


where each term in these equations is taken to be evalu¬ 
ated on the boundary. 


Taking = 0, we can then see that Eqs. (|A5|), (|A7|), 


(A8), (A9), and (AlO) allow the geodesic equation to be 


written as 

^ ^ - 2<1><1>, 


X^u = 


-( 4 ) 

Ht,x 7(3) 

- hlA 




tc=X 


+ 0(e6). 

(All) 


This is identical to Eq. (58). 


In this subsection we will follow the approach used 
by Chandrasekhar [STl l32] . If the 4-velocity is given by 
Eq. (37), then the components of the energy-momentum 


(Al) tensor are given by 


nab 


{p + pU+p)u‘^r+pg‘^\ 


(Bl) 


such that 

T“ = p{l+v^ + 11 + 2$) + 0(e®) 

T*^ = p ^1 + + n + 2$ + + O(e^), 


and 


T^^>' = pl^l+v^ +U + 2^ + ^j v>^v’' 

+ (1 - 2^)p5^'^ + 0(e®) . 

Let us now define a = p{l + + 2i> + 11 + p, and the 

total time derivative 


d d _ 


(B2) 


To derive the form of the conserved post-Newtonian mass 
let us consider 


r*% = cr,t + + p^,t - P,t ■ 


(B3) 


Using the continuity equations, and Eq. (B2), the last 
two terms in this equation can be written as 


dp f 

'^■•-’’■•=1’^-^-'’ A"-”-' 

^\dt^ 2 J dt 


(B4) 


Eq. (B2) can then be re-written as 


n=(:^ + v.u)u + q-($-> 


d 


dt 


d 


dt 


dp 

dt 


(B5) 

We can now use the continuity equations to relate the 
time dependence of the post-Newtonian energy density 
to the pressure: 

dn p dp 


p—— = - —- = —pV • V . 

^ dt p dt 


Thus, Eq. (B5) simplifies to 




(B6) 


(B7) 


Hence, at O(e^), we can use the conservation of energy- 
momentum to identify the following conserved post- 
Newtonian mass, Mpjsj'. 


MpN = 




+ 3$ 1 dV 


= ApA + ‘^{P^) ■ 


(B8) 
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